double a0 = 1e-6;
double b0 = 15e-6;

double a1 = 15e-6;
double b1 = 1500e-6;

double a2 = 15e-6;
double b2 = 1500e-6;

double ft(double *x, double *par) {

  double t = x[0];
  double tau_H = par[0];
  double tau_G = par[1];
  double N_H = par[2];
  double N_G = par[3];

  double tau_gen=par[5];
  double N_gen=par[6];

  double binW=par[7];

  double pH = N_H * TMath::Exp(-t / tau_H) / ((tau_H * (TMath::Exp(-a2 / tau_H) - TMath::Exp(-b2 / tau_H))));
  double pG = N_G * TMath::Exp(-t / tau_G) / ((tau_G * (TMath::Exp(-a1 / tau_G) - TMath::Exp(-b1 / tau_G))));
  double pGen = N_gen * TMath::Exp(-t / tau_gen) / ((tau_gen * (TMath::Exp(-a0 / tau_gen) - TMath::Exp(-b0 / tau_gen))));

  return binW*(pH + pG - pGen) + par[4];
}

int* h2dGetXY(int k){
    int *res=new int[2];
    res[1]=5-(k-1)/5;
    res[0]= ((k%5)==0)?5:(k%5);
    return res;
}

int draw() {

  gStyle->SetPaintTextFormat(".1f");
  TF1 *func_pol = new TF1("func_pol", ft, 1, 1.5e-3, 8);
  func_pol->SetNpx(1e6);
  func_pol->SetLineColor(2);
  func_pol->SetParNames("#tau_{LS}", "#tau_{Gd}", "N_{LS}", "N_{Gd}", "p0", "#tau_{Gen}", "N_{Gen}", "binW");
  func_pol->SetLineWidth(2);

  func_pol->SetParLimits(2,0,1e10);
  func_pol->SetParLimits(3,0,1e10);
  func_pol->SetParLimits(6,0,1e10);

  func_pol->FixParameter(0,210e-6);
  func_pol->FixParameter(1,30e-6);
  //func_pol->FixParameter(5,3e-6);
  func_pol->SetParLimits(5,1e-6,5e-6);

  double EdgesX[6]={0,1e6,2e6,2.5e6,3.25e6,4e6};
  double EdgesY[6]={-2000,-1400,-500,500,1400,2000};

  //TString top_dir = "/dybfs2/users/lijj/DATA/Accsub/IBD/IBD_version0.6/phi800mm_grid_time_neutron/Output/";
  TString top_dir = "phi800mm_grid_time_neutron/Output/";

  bool flag_Gd=false;

  int nAD[3] = { 2, 2, 4 };

  //fitted constant
  TH2D *h[3][4];

  TH2D *h_tau_Gd[3][4];
  TH2D *h_tau_LS[3][4];

  TH2D *h_N_Gd[3][4];
  TH2D *h_N_LS[3][4];

  for (int i = 0; i < 1; i++) {
    TFile *f = new TFile(top_dir + TString::Format("EH%d/EH%d_tot.root", i + 1, i + 1), "READ");
    for (int j = 0; j < nAD[i]/2; j++) {
      h[i][j]=new TH2D(Form("h_%d_%d",i+1,j+1),"",5,EdgesX,5,EdgesY);
      h_tau_Gd[i][j]=new TH2D(Form("h_tau_Gd_%d_%d",i+1,j+1),"",5,EdgesX,5,EdgesY);
      h_tau_LS[i][j]=new TH2D(Form("h_tau_LS_%d_%d",i+1,j+1),"",5,EdgesX,5,EdgesY);
      h_N_Gd[i][j]=new TH2D(Form("h_N_Gd_%d_%d",i+1,j+1),"",5,EdgesX,5,EdgesY);
      h_N_LS[i][j]=new TH2D(Form("h_N_LS_%d_%d",i+1,j+1),"",5,EdgesX,5,EdgesY);
      TH1D *h_voxel[25];

      for (int k = 0; k < 25; k++) {
        h_voxel[k] = (TH1D *)f->Get(Form("h1dDt_%d_grid_%02d", j + 1, k + 1));

        double binW=1e-6;
        if((k+1)==6 || (k+1)==7 || (k+1) ==8 || (k+1)==11 || (k+1)==12 || (k+1)==13 || (k+1)==16 || (k+1)==17 || (k+1)==18) {
            //h_voxel[k]->Rebin(10);
            flag_Gd=true;
            func_pol->ReleaseParameter(6);
            func_pol->SetParLimits(6,0,1e8);
            func_pol->SetParameters(210e-6, 30e-6, 10, h_voxel[k]->GetBinContent(2), 1, 3e-6, h_voxel[k]->GetBinContent(2),binW);
        }
        else {
            h_voxel[k]->Rebin(10);
            binW*=10;
            func_pol->SetParameters(210e-6, 30e-6, h_voxel[k]->GetBinContent(2), 10, 1, 3e-6, 10, binW);
            func_pol->FixParameter(6,0);
            //func_pol->ReleaseParameter(0);
            //func_pol->ReleaseParameter(2);
            //func_pol->FixParameter(1,1);
            //func_pol->FixParameter(3,0);
        }

        TFitResultPtr fptr=h_voxel[k] ->Fit("func_pol","SQE","",flag_Gd?1e-6:12e-6,1.5e-3);
        new TCanvas(Form("voxel_%d",k+1),Form("voxel_%d",k+1));
        h_voxel[k] -> Draw();
        cout<<Form("voxel %d, C/Ndf= %.2f",k+1,fptr->Chi2()/fptr->Ndf())<<endl;

        int iX=h2dGetXY(k+1)[0];
        int iY=h2dGetXY(k+1)[1];

        double tau_LS=func_pol->GetParameter(0)*1e6;
        double tau_Gd=func_pol->GetParameter(1)*1e6;
        double N_LS=func_pol->GetParameter(2);
        double N_Gd=func_pol->GetParameter(3);
        double constant=func_pol->GetParameter(4);

        double tau_LS_err=func_pol->GetParError(0)*1e6;
        double tau_Gd_err=func_pol->GetParError(1)*1e6;
        double N_LS_err=func_pol->GetParError(2);
        double N_Gd_err=func_pol->GetParError(3);
        double constant_err=func_pol->GetParError(4);

        if(N_Gd<1 || N_Gd<N_Gd_err) N_Gd=0;
        if(N_LS<1 || N_LS<N_LS_err) N_LS=0;

        h[i][j]->SetBinContent(iX,iY,constant);
        h[i][j]->SetBinError(iX,iY,constant_err);

        h_tau_Gd[i][j]->SetBinContent(iX,iY,tau_Gd);
        h_tau_Gd[i][j]->SetBinError(iX,iY,tau_Gd_err);

        h_tau_LS[i][j]->SetBinContent(iX,iY,tau_LS);
        h_tau_LS[i][j]->SetBinError(iX,iY,tau_LS_err);

        h_N_Gd[i][j]->SetBinContent(iX,iY,N_Gd);
        h_N_Gd[i][j]->SetBinError(iX,iY,N_Gd_err);

        h_N_LS[i][j]->SetBinContent(iX,iY,N_LS);
        h_N_LS[i][j]->SetBinError(iX,iY,N_LS_err);

        flag_Gd=false;
      }

      new TCanvas(Form("EH%dAD%d-Constant",i+1,j+1),Form("EH%dAD%d-Constant",i+1,j+1));
      h[i][j]->Draw("COLZ+TEXT+E");
      h[i][j]->SetMarkerSize(2);
      gPad->SaveAs(Form("Figures/EH%dAD%d-Constant.pdf",i+1,j+1));

      //new TCanvas(Form("EH%dAD%d-Tau-Gd",i+1,j+1),Form("EH%dAD%d-Tau-Gd",i+1,j+1));
      //h_tau_Gd[i][j]->Draw("COLZ+TEXT+E");
      //h_tau_Gd[i][j]->SetMarkerSize(2);

      //new TCanvas(Form("EH%dAD%dTau-LS",i+1,j+1),Form("EH%dAD%d-Tau-LS",i+1,j+1));
      //h_tau_LS[i][j]->Draw("COLZ+TEXT+E");
      //h_tau_LS[i][j]->SetMarkerSize(2);

      new TCanvas(Form("EH%dAD%d-N-Gd",i+1,j+1),Form("EH%dAD%d-N-Gd",i+1,j+1));
      h_N_Gd[i][j]->Draw("COLZ+TEXT+E");
      h_N_Gd[i][j]->SetMarkerSize(2);
      gPad->SaveAs(Form("Figures/EH%dAD%d-N-Gd.pdf",i+1,j+1));

      new TCanvas(Form("EH%dAD%d-N-LS",i+1,j+1),Form("EH%dAD%d-N-LS",i+1,j+1));
      h_N_LS[i][j]->Draw("COLZ+TEXT+E");
      h_N_LS[i][j]->SetMarkerSize(2);
      gPad->SaveAs(Form("Figures/EH%dAD%d-N-LS.pdf",i+1,j+1));
    }
  }

  return 0;
}
